function idat = fsb_detrend_scan(idat,sandbox,thresh,detopt)

% FSB : remove a linear drift from scanner signal
%
% EXAMPLE:
% idat = fsb_detrend_scan(idat,sandbox,0,0)
%
% INPUT:
% idat: 4-D image data
% sandbox: sandbox experiment struct
% thres: threshold setting
% detopt: 0 linear, 1 nonlinear
%
% OUTPUT:
% idat: detrended 4-D image data
%
% CALLED BY:
% FSB.m
%
% NOTES:
% Contrast problem in detrending fixed now.
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
%$Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  explicit = 1; % Give graphics output

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Check if correct dataset size
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

d = size(idat);
if size(d)<4;
    disp('No 4D dataset, aborting');
    return
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Look for threshold input, resort to default if missing
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if nargin<3;
    thresh = 1000;
    detopt = 0;
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculate idat change values before detrending
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

idat2 = single(idat);
idat2(idat2<thresh)=NaN;
avgint = fsb_meancalc(idat2,4);
sub_disp_image_mean(idat2,'all','Before processing scans');


if ~exist('sandbox','var') && ~isempty(sandbox.intrial) && size(sandbox.intrial,2)>3;

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Run function without a sandbox input
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    disp('detrending ...')
    h = waitbar(0,'Detrending...');
    for x = 1:d(1);
        waitbar(x/d(1));
        for y = 1:d(2);
            for z = 1:d(3);
                te = double(squeeze(idat(x,y,z,:)));
                meante = mean(te);
                idat(x,y,z,:)=int16(detrend(te)+meante);
            end
        end
    end
    close(h);
    disp('...done');
    hold on;
    avgint2 = fsb_meancalc(idat2,4);

else

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Run function with sandbox input
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    g = waitbar(0);
    if size(sandbox.intrial,2)<4
        sandbox.intrial(:,4)=1;
    end
    for scannum = 1:max(sandbox.intrial(:,4));

        detstring = ['detrending scan ' num2str(scannum) ' of ' num2str(max(sandbox.intrial(:,4)))];
        disp(detstring)
        waitbar(scannum/max(sandbox.intrial(:,4)),g,detstring);

        scanindex = find(sandbox.intrial(:,4) == scannum);
  
        idat_seg = single(idat(:,:,:,scanindex));
        intrial_seg = sandbox.intrial(scanindex,:);

        dim = size(idat_seg);
        idat3 = shiftdim(idat_seg,3);
        idat3 = reshape(idat3,size(idat3,1),[]);

        if detopt == 0

            meante = nanmean(idat_seg,4);
            meante = meante(:);
            meante = repmat(meante,[1 length(scanindex)]);
            detdata = detrend(idat3);
            det = meante'+detdata;

        else

            idat3(:,find(intrial_seg(:,1)==0)) = NaN;
            det = fsb_detrend_nan(idat3,3,0);

        end

        idat_seg = reshape(det,[dim(4),dim(1),dim(2),dim(3)]);
        idat_seg = shiftdim(idat_seg,1);
        idat_seg = int16(idat_seg);

        disp(['...done with scan ' num2str(scannum)]);
        idat_seg_red = idat_seg;
        idat_seg_red(idat_seg_red<1000)=NaN;
        sub_disp_image_mean(idat_seg_red,scannum,'After detrending');
        idat(:,:,:,scanindex) = idat_seg;

    end

    avgint2 = fsb_meancalc(idat,4);
    close(g);
end

idat2 = idat;
idat2(idat2<thresh)=NaN;
sub_disp_image_mean(idat2,'of all','After detrending of all scans');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Create plots
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if explicit ==1;
    
    figure(501);
    set(gcf,'Name',' Plot of mean volume intensity before (blue) and after (red) Detrending');

    plot (avgint); hold on;
    plot (avgint2,'r');
    xlabel('Volume in trial');
    ylabel('Mean volume intensity');
    legend('Before','After')

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Subfunction to calculate mean volume change over time
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    function sigchange = sub_disp_image_mean(idat,scannum,scanstr)

        avgint3 = fsb_meancalc(idat,4);
        avgintn = smooth(avgint3);
        figure(701);
        set(gcf,'Name','Smoothed volume timecourse');cla;
        plot(avgintn);
        firstint = mean(avgintn(1:round(length(avgintn)/4)));
        lastint = mean(avgintn(round(length(avgintn)/4):end));
        sigchange= (lastint-firstint)/firstint*100;
        clear firstidat lastidat;

    end


end
